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PROJECT SUMMARY 


This research was aimed at developing a novel concept for optimally designing output feedback 
controllers for plants whose dynamics exhibit gross changes over their operating regimes. This was, 
essentially, to formulate the design problem in such a way that the implemented feedback gains 
vary as the output of a dynamical system whose independent variable is a scalar parameterization 
of the plant operating point. 

The results of this effort include derivation of necessary conditions for optimality for the general 
problem formulation, and for several simplified cases. The question of existence of a solution to 
the design problem was also examined, and it was shown that the class of gain variation schemes 
developed in this effort are capable of achieving gain variation histories which are arbitrarily close 
to the unconstrained gain solution for each point in the plant operating range. The theory was 
implemented in a feedback design algorithm, which was exercised in a numerical example. 

The results of the research undertaken under this contract are applicable to the design of 
practical high-performance feedback controllers for plants whose dynamics vary significantly during 
operation. Many aerospace systems fall into this category. 
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1. INTRODUCTION 


Linear control theory has provided many techniques for the design of feedback controllers for 
plants with linearizable dynamics. In most of these methods, the feedback synthesis is performed 
for a plant model obtained by linearizing about some “representative” setpoint. Modelling errors 
and minor parameter variations are accommodated by taking measures to ensure that the control 
design is “robust.” A robust controller, for our purposes, is one which performs adequately when 
implemented with plant dynamics which differ somewhat from the nominal design model, in the 
sense of their being a perturbation of the nominal plant. Robust design methods have been an 
active area of research in recent years [1 >4], and have been brought to a fairly high level of 
maturity. 

Many dynamical systems encountered in practice, however, exhibit gross, structured changes 
in their dynamics as they move about in their operating regimes. For these systems, it is difficult 
to design feedback controllers with globally satisfactory performance using design theory intended 
to accommodate only local plant perturbations. This problem motivated work which led to a num- 
ber of approaches for using information about global plant parameter variations in designing the 
feedback gains. References [5 - 7] develop a design procedure in which a vector-valued cost function 
reflecting the plant dynamics at a number of operating points is Pareto-optimized. Implementation 
of the method leads to a nonlinear programming algorithm. In [8, 9] a design approach is described 
in which loci of permissible gain values are established for each of several operating points. Per- 
missibility, here, means that the gains for each model place its closed-loop poles in some specified 
region of the complex plane. The feedback design is performed through examining tradeoffs over 
the intersection of these loci. References [10, 11] describe quadratic optimization procedures for 
stabilizing the members of a discrete set of linear plants by using output feedback gains which 


minimize a scalar cost function reflecting the performance of all of the plants in the set. Finally, 
in [12 - 15], the problem of stabilizing a set of plants with a single compensator is examined in 
an algebraic context. A number of elegant results are derived, including global parameterization 
of stable compensators which stabilize a particular plant and, conversely, plants stabilizable by a 
given stable compensator [12, 13]. This approach to the global feedback design problem does not, 
however, appear to be developed to quite the level of practical applicability seen in the methods 
based on quadratic optimization [14, 15]. 

There is a serious practical difficulty with all of the work described above. The compromises 
required to make a fixed set of feedback gains perform adequately over a system’s entire operating 
range can result in less satisfactory closed-loop performance at any single operating point than that 
attainable using a feedback scheme which can vary to accommodate changes in the plant. These 
latter feedback schemes fall into two broad categories. The first includes self-tuning regulators 
[16, 17] and controllers incorporating “online redesign” logic [18]. The second category consists of 
feedback structures whose gain elements are designed to have a particular functional dependence 
on parameters coordinatizing the instantaneous system operating point. The most typical of these 
latter approaches is gain scheduling. 

In situations where the plant dynamics are relatively well-known, gain scheduling and similar 
gain variation strategies have an important conceptual advantage over those employing online tun- 
ing. The fact that the actual controller design is performed offline permits the designer as many 
testing and redesign iterations as are necessary to achieve the desired performance characteristics 
throughout the plant’s operating range. Another consideration, which becomes important when 
there is a premium on the controller implementation’s computational overhead, is that the structure 
of the gain’s dependence on the operating point parameterization is chosen by the designer. This 
permits further flexibility in deciding the tradeoff between the controller’s performance and the 
complexity of its implementation. These considerations make this class of gain variation scheme es- 
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pecially well-suited for aerospace applications and it is, in fact, the one most commonly encountered 
in aerospace industrial practice. 

Neither of the two gain scheduling approaches common in industrial practice exploit these 
advantages efficiently. The first approach consists simply of defining a grid over the domain of the 
plant parameterization and designing a set of feedback gains for each grid point. In implementation, 
the gains are obtained by interpolating to the current operating point from these grid points. This 
brute force approach can be very expensive in terms of the control computer’s memory requirements. 
The other approach is to implement the gain as some simple curve fit, often straight line, to 
gains calculated at selected parameterized operating points. This cures the problem of excessive 
computer storage requirements, but introduces grave difficulties of its own. The feedback gains thus 
implemented are only approximations of the designed gain values, and there is no a priori guarantee 
that the approximate gains will perform well. The practical result of this is that designing gain- 
scheduled feedback in this piecemeal fashion characteristically involves a great deal of uncertainty 
and “cut-and-try.” 

Recently, the quadratic optimization procedure developed in [11] has been extended to provide 
a much more orderly and rigorous procedure for designing scheduled-gain feedback [19, 20], Essen- 
tially, the extension consists of redefining the plant models reflected in the problem cost function 
in such a way as to embed the gain schedule structure in their input and/or output matrices. This 
method can accommodate any gain schedule structure expressible as a polynomial in the param- 
eters used to coordinatize the plant operating points. Reference [21] provides a full derivation of 
the procedure, and describes a nontrivial application to a self-repairing flight control problem. 

The principal - and significant - advantage of this approach over defining gain schedules though 
the use of approximate curve fitting methods is that the scheduled gain is “exact” at each designed- 
for operating point, in the sense that the optimization is based on the scheduled gain’s actual effect 
on the closed-loop dynamics of each plant model appearing in the design. Because of this, at 
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least at the operating points used in designing the controller, the effect of the scheduled feedback 
in implementation will be exactly as designed for, modulo effects due to modelling errors. The 
method’s difficulties, on the other hand, are twofold. First, as in [10, 11], a discrete set of linear 
plant models is used to represent the plant’s globed parametric variations. The choice of these 
models is, essentially, a matter of heuristics and experienced judgment. Second, although the 
method can be used to design polynomial gain schedules of arbitrary finite order, the number 
of free parameters to optimize increases at least linearly with polynomial order unless further 
restrictions, again based on heuristics and judgment, are imposed on the form of the schedule. 

This report develops an quadratic optimization-based approach for designing feedback gains 
which vary with system operating condition, applicable when the plant’s location in its operating 
regime can be parameterized by a scalar function. The gain matrix is designed to be the output of 
a dynamical system having the plant operating point as its independent variable. This is a class of 
gain variation constraints which contains gain scheduling as a subset. In gain schedules, the gain 
is a function only of the independent variable; i.e., the system operating condition. The approach 
developed in this report can be used to design gain variation schemes in which that restriction is 
relaxed, to permit functional dependence on the instantaneous gain “state.” Because of this, the 
gains can be designed to vary over the plant’s operating regime in ways which would require infinite 
series representations, if implemented as polynomial gain schedules. On the other hand, the order 
and structure of the gain variation dynamics can be chosen by the designer in a tradeoff between 
computational overhead in the controller implementation, and the plant’s requirement for complex 
gain variations. An additional positive feature of the theory developed in this report is that, as 
a direct consequence of the problem formulation, the design reflects the plant dynamics as they 
vary continuously across the operating regime, rather than only at an arbitrarily chosen collection 
of setpoints. In the sequel, we will refer to this class of gain variation schemes as PDGP schemes, 
short for Parameter-Dynamic Gain Propagation. 
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Section 2 first formulates the PDGP feedback optimization problem with only minimal as- 
sumptions made on the form of the gain propagation dynamics, then proposes and examines several 
candidate PDGP constraint expressions having features which make them attractive for practical 
implementation. The Section also briefly examines the question of existence of a solution to the 
general optimization problem. Section 3 derives necessary conditions for optimality for one of the 
PDGP constraint structures examined in Section 2, and uses these to develop a numerical design 
algorithm, which is exercised in a numerical example. Conclusions are presented in Section 4. 

A significant byproduct of the numerical algorithm development has been the formulation of 
a novel, efficient iterative procedure for solving discrete Lyapunov equations. A key feature of the 
theory leading to this method is a least-squares-optimal approximation to certain symmetric sums 
of Kronecker products. This development is documented in Appendix A. The Lyapunov equation 
algorithm and its derivation are provided in Appendix B. 
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2. PDGP PROBLEM FORMULATION 


This Section formulates the general PDGP optimization problem on which the numerical design 
algorithm of Section 3 is based. We consider a plant with input u € Z m and output y € Z L described 
by 


d_ 

dt 


x 


f(x,u,a) + w(t,a) 


( 1 ) 


y = y(x,a) + v(t,a) (2) 

where x € Z n is the plant state, and w and u are white noise processes. The “plant parameter” a 
gives the location of the plant in its operating regime. It reflects quantities which vary slowly enough 
during plant operation to be modeled as constant in the feedback design, but whose values affect 
the plant’s dynamic response. Examples of quantities which are often chosen as plant parameters 
include aging in process control, and Mach number in flight control problems. For this report, it 
is assumed that a varies on the domain a 0 < a < a/. Assuming that the plant dynamics are 
linearizable throughout the domain of a, assign a locus of setpoints x and u as a function of a; that 
is, {x(a), u(o)} : Z —* Z n x R m . The plant perturbation dynamics are then described in discrete 
time by 


x(k -f l,a) = A(ot)x(Jk,a) + B(a)u(k,a) + tu(fc,a) 


( 3 ) 


y(fc,ot) = C(a)x(fc,a) + v(k,a ) 


( 4 ) 
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where k is the time sample index. In (3, 4), the variables w and v are discretizations of w and v 
with 


E{w(k, ct)v> T (a, a)} = W (a)$fc, (5) 

E{v{k,a)v T (a t a)} = V(a)Sk, (6) 

£?{u;(*,a)v T (s,a)} = 0 (7) 

The plant is to be stabilized for all a € [a 0 ,o/] by an output feedback control law 

u{k, a) = —G(a)y(k, a) (8) 

in which the feedback gain varies with a according to 

« = r(0,a) (9) 

<2 = 11 ( 0 , a) (10) 


where the “dot” notation on 0 denotes differentiation with respect to ot rather than time. Free 
parameters are 0(d), the value of 0 E Z q at some fixed d in the domain of a. The PDGP feedback 
design problem, then, consists of choosing 6 and the functional form of T and II, then adjusting 
the value of 0(d) to achieve the desired performance. The potential complexity of this problem 
motivates the use of optimization methods in the design process. 

A cost function commonly used in discrete-time optimal output feedback design at a given 
operating point [11] is 
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c(a) = lira 
v N-+00 


2(tf+l) 


N 

^2 E{x T (k + 1, a)Q(a)x(k + 1, a) + u T (k, a)R(a)u(fc, a)} 


fc =0 


( 11 ) 


where Q(a) > 0, Jt(a) > 0. Reference [11] develops a number of sufficiency conditions for the 
existence of a G which minimizes (11) at a given operating point. Thinking of c(a ) as the “a 
- instantaneous” system performance, it is reasonable to define a global cost over the operating 
regime by 


J = f c(a)dot (12) 

J a« 

In the interest of clarity, notation indicating dependence on ot will be suppressed in the sequel. 
The problem of minimizing (12) subject to (3, 4) and (9, 10) is more conveniently posed as that of 
minimizing a Lagrangian function. It is well known that, when the closed-loop dynamics defined 
by (3, 4) and (8) are asymptotically stable, c can be expressed as 

c = tr{KW} + tr{G T {B T KB + R)GV} (13) 

where K satisfies the constraint 


$ {K, G) = <f> T {G)K<j>(G) - K + Q + C T G T RGC = 0 


(14) 


<I>(G) -A- BGC 


(15) 


The constraints (9, 10) and (14) are adjoined to J, using the expression (13), to form the Lagrangian 


£= r'c + <r{J(K,G)Aj}+lr{[G-n(«)lAS} + Anr(*)-«l*» ( 16 ) 

where Ajf> Aq and A $ are Lagrange multipliers associated with the constraints on K, G and 9. 
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Denoting the integrand of (16) as X, the Euler-Lagrange equations for a stationary point of L 


are 


dX/dK = 0 d)//dA K = 0 


dX/dG = 0 dX/dA a = 0 


A t = -dX/d0 O = T{0,a) 


Equations (17.a), expanded, are discrete Lyapunov equations: 


dX/dK = k4 T - A K + W + BGVG t B t = 0 


dU/dA K = S(K,G) = 0 


Equations (17.b), expanded, are 


dX/dG = kGA K - B t KAK k C t + A 0 = 0 


Z=B T KB+R 


A K = CA K C T + V 


dX/dtia = G - 11(0, a) = 0 


(17.a) 

(17.b) 

(He) 

(18) 

(19) 

(20.a) 

(20.b) 

(20.c) 

( 21 ) 
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and the equation for A* from (17.c) is 


A, = ^H n (*,a)AS}-Aj- r( *,a)l (22) 

There remains the problem of optimally assigning boundary conditions on 9 and There is no 
terminal cost in this problem, so A# (a/) is chosen as 

A,(a/) = 0 (23) 

The effect of the choice of 9(a 0 ) on the cost is obtained by integrating X T 0 by parts in (16), which 
leads to the necessary condition 


d£/d9{a 0 ) = \ 9 {a o ) = 0 (24) 

In summary, then, once a structure for T and II is chosen, the PDGP optimization problem 
reduces to a matter of determining the optimal initial condition for 6. Since only the initial 
conditions are free in this problem, the key to effectively employing the above theory in designing 
a feedback gain variation scheme is to hold elements of 9 constant, thus assigning them the role 
of static free parameters in the optimization. For example, a typical single-variable gain schedule 
takes the form 


G(a) = N + (a - a 0 )M 
This can be reexpressed as a PDGP scheme: 


(25) 


G = M G(a„) = N 

With the above notation, one could choose 9 as 


(26) 
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9 = vee[G : M) 


(27) 


where vec(.) is the column-stacking operator, so that 


T = vec[M : 0] 

(28) 

n = mat{[I m i : Q]9} 

(29) 


where mat(.) denotes column unstacking and J m i in an identity matrix of dimension ml. Given the 
simplicity of (28, 29) however, it is easier to separate G and Af , and assign Lagrange multipliers 
A g and Am- The Lagrangian for this problem, then, is 

Z = [*' c + tr{S(K,G)hl ) + tr{[M - G) a£) + tr{-MA^}da (30) 

Jo t 0 

and the Euler-Lagrange equations are (18, 19) and 


G = M 

(31) 

M = 0 

(32) 

-dX/dG = Ag = B t KAK k C t - KG\ k 

(33) 

-dX/dM = A m = -Ag 

(34) 


with boundary conditions 
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Ao(a 0 ) = A a (or/) = 0 


(35) 


Awfoio) = Am (a/) = 0 (36) 

We now turn our attention to PDGP structures in which the gain matrix varies as a linear 
combination of states satisfying a system of linear homogeneous ordinary differential equations 
(ODEs). Two of the simplest T structures suitable for implementation in a flight control computer 
are 


9 = N9M (37) 

9 =N9 + 9M (38) 

where 0 6 R p * q , N € Z pXp and M G Z qXq . These are both special cases of 

vec 9 = F vec 9 (39) 

Their implementations! advantages over (39) lie in storage and computational considerations. In 
(39), F G Z pqXpq and p 2 g 3 multiplications are required to calculate 9. In (37) and (38), N and M 
occupy p 1 + q 1 locations, and the derivative calculation requires gp 5 + pg 1 multiplications. This 
difference becomes important for large-order PDGS schemes. A suitable choice of II structure for 
either (37) or (38) is 


G = S9D + P 


( 40 ) 


The elements of S, D and P can be either free parameters in the optimization, or fixed. 
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The chief characteristic of linear PDGP structures such as those defined in (37 • 40) is their 
analytical simplicity. This by no means implies that these structures are limited in the level of 
closed-loop performance they can attain, however. The following lemma quantifies this issue. 
Lemma 1 : Assume that A(or), B(a), C(a), Q(a), R(ct), W (at) and V(a) have continuous k tK 
derivatives on the interval [ot 0i aj) for k > 0. Further, assume that the system (3, 4) is stabilizable 
by output feedback and that for each a €E [at e ,a/] there exists a 

G* (a) = arg min {c(a)> (41) 

Qt${a) 

for c(a) defined in (11), and where p(a) is the set of asymptotically stabilizing gains for the plant 
at the operating point parameterized by the current value of a. Then, for 9 € Z q given by (39) 
with q > 2 mlr for integer r > 0, and G given by 

vec G = S vec 9 + P (42) 

there exist F , S, P and 0(a o ) such that 

max ||<T (a) - G(a)|| < b/(r + l) fc+1 (43) 

ac[a.,a/] 

where b > 0 is a constant. 

Proof : See Appendix C. 

The basic implication of this lemma is that by specifying gain dynamics of sufficiently large 
dimension in a linear PDGP problem, one can construct a gain variation which is arbitrarily close 
to the unconstrained variation of G* over the operating regime. We note that the bound (43) was 
obtained by constructing an individual Fourier series expansion for each element of G and is, thus, 
extremely conservative. 
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In the next Section, the PDGP optimization problem will be solved for a slight simplification 
of the structure (38, 40). One of the attractive features of (38) is that it leads to a very simple 
transition expression: 


$ (a + A a) = e NAa 0(a)e MAa (44) 

This can be easily exploited to permit use of the theory developed thus far in designing controllers 
for systems whose operating regimes are best coordinatized by a vector-valued parameter. For 
example, consider the case of an operating regime parameterized by two variables, say, pi and pj. 
Now, define a locus of p = (pi(a), p 2 (a)) in the operating regime forming a path, with arc length 
a as the independent variable. This situation is illustrated in Figure 1, for p chosen to be a spiral. 
In order to move from point A to point B, one simply integrates (38) between ot(A) and ct(B). To 
move to point C, rather than laboriously propagate (38) around three revolutions of the spiral, it 
would be more economical to have e NAa and e MAa for Aat roughly corresponding to one full turn 
stored in the controller implementation. This would permit using (44) to “cut across” turns of the 
spiral to the turn containing C. The final correction in reaching C would then be performed, again 
using (38). 
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3. A LINEAR PDGP SCHEME 


In this Section, necessary conditions for optimality in a simple linear PDGP structure are 
derived. These, in turn, are used in a computational algorithm for calculating locally optimal gain 
parameters. The algorithm is exercised in a numerical example. The class of gain propagation 
schemes considered take the form (38, 40) with S and D in (40) restricted to identity matrices, 
and P = 0; in other words, 

G = NG + GM G(a 0 ) = G 0 (45) 

N = 0 (46) 

M = 0 (47) 

The Lagrangian for this problem is 


L= [*' e + triSiK'G^fo + trUNG + GM- 6} K%) 

Jet* 

+ tr{-MAj[ f } + tr{-//Aj))<ia (48) 


and the Euler-Lagrange equations are (14) and (18), (45 - 47) and 


A 0 = — d)i/dG = - N t A g ~ A g M t - KG\ k + B t KAA k C t (49) 


A n = — dX /dN = — A qG^ 


(50) 
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ku = -dX/dM = -G T A a 


(51) 


The boundary conditions to be satisfied are 


Aa(a 0 ) = Aa(otf) = 0 

(52) 

Aft(<*o) = A *r(<*/) = 0 

(53) 

Am(«o) = A m(<»/) = 0 

(54) 


Solution of the necessary conditions for optimality in this problem consists, essentially, of 
determining G 0 , N and M such that A a. Am and As have null initial conditions. The simplest 
general solution approach is to exploit the fact, from (24), that 


dZ/dG 0 = A 0 {ct 0 ) 

(55) 

dC/dN = A N (a 0 ) 

(56) 

dZ/dM = Am(«o) 

(57) 


Equations (55 - 57) suggest converting the problem of satisfying the costate boundary conditions to 
one of locally minimizing the Lagrangian using a descent algorithm. For example, a simple steepest 
descent algorithm takes the form: 

0. Determine an initial guess for G°(ctf) t N° and M° . Set k — 0. 

1. Integrate (49-51) backwards from a/ to a 0 , to obtain A^(a 0 ), A%(a 0 ) and Aj^(a„). 
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2. Increment G 0 , N and M as 


G *+i = G k_ CfcA * (« 0 ) (58) 

N k+i = N k_ e fc A^(a 0 ) (59) 

M<* +1 > = Mfc - £A*,(a 0 ) (60) 

. where £jt € (0, 1] is chosen to ensure a satisfactory cost decrease for the iteration. One standard 
criterion, given in [22], for determining (k takes the following form for this problem: 

P$k<Pk <J k - J k+1 < (1 - p)tk<Pk (61) 

. where p is a fixed parameter in the interval (0,0.5], J k is the value of J from (12) for the k th 
iteration’s parameter values, and 

P* = IIAJMII* + »AjMa.)||* + ||Aj,( a .)|| J (62) 

3. Set k = A: + 1 and go to 1. 

In order to put the above theory in perspective, we now consider a feedback design problem 
for a simple linear system with variable dynamics, given by 

* = A c + u + «/ (63) 

y = * (64) 

where A e is 
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A — \ 0 1 

l"" 3 («) ~ 01 


( 65 ) 


and a/* varies as 

o/*(a) = 0001 + ( 10 - «X)l)a (66) 

for 0 < ot < 1. The system is sampled at 10 Hz, and it is assumed that the discretized process 
noise covariances are W = I and V = 0. The penalty weights are Q = jR = /. 

Feedback gains were calculated using the PDGP formulation of this Section and, for compar- 
ison, the design approaches described in [11] and [19 - 21]. In [11, 19 - 21], the cost function is 
defined for a discrete collection of plant models, which are chosen by the designer to represent plant . 
variation over the domain of a . The cost takes the form 

•/=£>(“»> (« ? ) 

j=i 

where n p is the number of models in the collection, and c(oty) is the quadratic regulation cost 
defined in (11). For the example problem, 11 plant models were chosen, spaced over 0 < a < 1 at 
intervals of Aat = 0.1. 

The necessary conditions for optimality consist of satisfying (18, 19) at each of the otj and 

KjG{ k K )j ~ Bj KjAjMjCf = 0 (68) 

j=i 

where the j subscripts denote evaluation at otj. Equation (68) can be viewed as a discrete-parameter 
analog of (20.a). In [11], G does not vary with a. References [19 - 21] develop an extension of the 
theory in [11] in which the gain schedule structure 

G[a) = G 0 + aG t (69) 
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is embedded into the optimization problem formulation. It should be noted that the theory in [11, 
19 - 21] also accommodates vector-valued a, in addition to the scalar case. 

Optimal feedback solutions for the fixed-gain, scheduled-gain and PDGP formulations are given 
in Table 1. In order to represent unconstrained variation of the optimal gains, “pointwise” optimal 
output feedback gains (n p = 1) were calculated for each of the plants in the fixed and scheduled- 
gain design model. The numerical calculations were considered converged when the 2-norm of the 
gradients Ao(0) or (68) dropped below one percent of the cost. In order to consistently compare the 
performance of the various feedback variation structures, the a-integral performances (12) of the 
fixed, scheduled and pointwise-optimal designs were calculated from their individual model costs, 
using Simpson’s rule. 

Unsuprisingly, the fixed-gain design gave the worst performance, 14.25, whereas the uncon- 
strained pointwise-optimal design returned a cost of 13.24. For this example, the PDGP and 
scheduled-gain designs returned 13.38 and 13.36, respectively. These values can be considered vir- 
tually the same, given the rather casual accuracy of the Simpson’s rule quadrature used to obtain 
the cost for the scheduled-gain design. The variation of the model cost c(a) for each of the gain 
solutions is displayed in Figure 2. The Bcheduled-gain design optimization apparently sacrificed the 
performance at or = 0 in order to achieve a good match to pointwise optimal performance across 
0.1 < a < 1. The PDGP scheme, on the other hand, resulted in less severely degraded performance 
at a = 0, but allows the degradation to persist until approximately a = 0.12; thereafter, it very 
slightly outperforms the scheduled gain. 

Figure 3 displays the gain variation histories for the elements of G. It would appear that 
the 1,1 element of G is bears primary responsibility for adding damping to the system as in 
(65) increases. Both the PDGP and scheduled-gain solutions closely approximate the trajectory 
of this element. The scheduled-gain and PDGP histories for the other elements, however, are 
significantly different, particularly for the off-diagonal terms. It must be remembered that the goal 
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TABLE 1. FEEDBACK GAIN VARIATION SOLUTIONS 


FIXED: 


1.9149 - 0 . 1468 ' 
- 0.0845 0.6476 


SCHEDULED: G = G a + aG x 


Go 


0.9263 0.0471 
0.1670 0.9882 


G i = 


1.7343 - 0.2920 
- 0.8645 - 0.7316 


PDGP 


AC* 

~P- = NG + GM 
da 


G(0) = 


1.0693 0.0057 
- 0.0076 0.8944 


0.4188 - 0.2062 

0.1956 3.5646 


M = 


0.5246 

- 0.0015 


0 . 3326 ' 
- 4.4935 


POINTWISE COST 







of the scheduled-goal of the scheduled-gain and PDGP methods considered here is not necessarily 
producing gains which provide a “best-fit” to the unconstrained variation of optimal G , but rather, 
to optimize the system performance over the domain of a a, subject to the gain schedule or PDGP 
dynamic constraints. On the other hand, as we know from Lemma 1, if a PDGP structure with 
more degrees of freedom had been chosen, the trajectory matches between the unconstrained and 
PDGP gains would have been closer. 
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4. CONCLUSIONS 


This report has documented a research effort which developed and examined a novel optimization- 
based approach to designing multivariable output feedback gains which vary to accommodate 
changes in the plant dynamics being regulated: Parameter-Dynamic Gain Propagation, or PDGP. 
The key feature of the approach is that the designed gains vary as the output of a dynamical system 
which, as its independent variable, uses a parameter which coordinatizes the plant dynamics in a 
manner analogous to the parameter used in traditional gain scheduling schemes. 

The general optimization problem was formulated, and necessary conditions for optimality 
were derived. Conditions for the existence of a solution to the problem were also obtained. Several 
simplified PDGP structures were examined, with special attention being paid to a restricted version 
of the general linear structure, for which a numerical optimization algorithm was derived and 
implemented. This latter was tested and demonstrated in a numerical example. 
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APPENDIX A 


The work described in this appendix was performed under NASA Contract NAS1-17493 and 
is included here for completeness. 
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APPENDIX A. 


LEAST SQUARES APPROXIMATION 
OF SYMMETRIC KRONECKER PRODUCT SUMS 

A situation which occurs with some frequency in optimization analyses involving matrix-valued 
quantities is that of having to invert positive definite matrices having the general form 

o v 

H = + [Cj ® Cj + Cf ® Cf] (A.l) 

<=i j= i 

where A\ > 0, JB t - > 0 for * = 1, and Cj , j = 1, ..., u are general square matrices. Numerical 

inversion of H is usually very cumbrous unless the terms in the Kronecker products are of relatively 
small dimension. Furthermore, one of the most common situations for which a requirement for 
H~ l exists is in Newton or quasi-Newton iterative schemes, where H is the Hessian of a scalar cost 
function which is to be minimized. In these cases, the exact value of H~ x is not actually needed, 
insofar as the inverse Hessian is only used to exploit “curvature” information and to provide scaling 
in modifying the gradient-based descent step in the problem’s free parameters. This motivates the 
notion of approximating H by a quantity H, having the structure 

H = S®P (A.2) 

so that its inverse is 


H~ l = S " 1 ® P~ l 


(A.3) 


When the Cj are zero, a straightforward approach to solving this approximation problem is to 
calculate S and P to minimize the Euclidean norm of H — S ® P, given by 
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J = tr{H T H - (S r ® P T ) H - H t (S ® P) + S T S ® P T P} 
The necessary conditions for a minimum in J are 


(A.4) 


dJ/dS = 0 dJ/dP = 0 


(A.5) 


or 


n 

dj/as = ||P|| J S - Aitr{BiP} = 0 (A.6) 

i=l 


dJ/dP = ||S|| a P - £ = 0 (A.7) 

»=i 

Equations (A.6) and (A.7) can be solved iteratively, using the following simple algorithm: 

0. ) Choose any pl°) > 0. Set k = 0. 

1. ) Set k = k + 1. 

2. ) From (A.7), evaluate 


S<*> = 


1 


• sl 


(A.8) 


3.) From (A.8), evaluate 


p(*) 




BitriAiS^} 


(A.9) 


4.) Go to 1. 

Numerical experience [21] indicates that this algorithm is quite robustly convergent, and a formal 
proof of algorithmic convergence is under development. It can easily be seen from (A.6) and (A.7) 


that 
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d t J/dvecSdvec T S = ||P|| 3 J n a Xn » (A.10) 

d 2 J/dvecPdvec T P = ||S|| 3 / n9xna (A.ll) 

so that (A. 8) and (A.9) are a quasi-Newton iteration. 

In order for the approximation (A.2) to be of practical value, it is important that S > 0 and 
P > 0, so that the inverse in (A.3) is positive definite. Furthermore, since the algorithm is only 
exercised through a finite number of iterations, it is desirable to know for what k one can guarantee 
that >0 and pW > 0. These issues are settled in the following lemma: 

Lemma2 : If H > 0 then S ^ > 0 and P^ > 0 for all A: > 0. 

Prop/ : See Appendix C. 

Unfortunately, a direct extension of the above theory to the case where the C 3 - are not zero 

does not generally result in H > 0. In order to treat this case, it is necessary to constrain S > 0 

and P > 0. This is done for Aj, B% and Cy € £ nXn by enforcing 


S = §S T P = PP T (A. 12) 

where S € £ nXn * and P € Jl nXn * , where n, and n p may be freely chosen. For this case, we seek 
a minimum of 


/ = tr{H T H - BS T ® PP t H - H t SS t ® PP T + SS T SS T ® PP T PP T > (A.13) 


The necessary conditions for a local minimum of J are 
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dJ/dS = 4 tr{(P T P)*}SS T S 

-i^tripTBiPyAiS 

» = 1 

-4^T tr{P T CjP} [Cj + Cj)S = 0 ( A.14) 

y=i 


d J/dP = 4tr{(S T 5)*}PP T P 

-<f< r (5 r AW 

»=l 

- 4 2 /r{5 r C J 5>(C/ + C,)P = 0 (A. 15) 

3-1 

Equations (A.14) and (A.15) could be incorporated into a steepest descent algorithm for performing 
the minimiz ation of J. This is not particularly attractive, however, because of the known slow 
convergence of steepest descent methods. It is preferable, instead, to seek a numerically inexpensive 
means of incorporating some, if not all, of the second derivative information into the iterations. 
The simplest is to use the increments 


6§ij = 


dJ/dSij 

d*J/dS £ 


SPa = 


dJ/dPij 

d'J/dPij* 


(A. 16) 


where (.),•/ is the ij tk element of the matrix. This is equivalent to using the inverse of the main 
diagonals of d 1 J/dvecSdvec T S and d 7 J/dvecPdvec T P in the step calculations. This approach has 
been used in other general optimization problems [23], and has provided a significant enhancement 
in algorithmic performance over steepest descent in those applications. The main diagonal elements 
of the Hessian matrix are given by 
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4tr{(P r P) J }[(5 r 5) iJ + (SS T ) ij + 5?] 


d*J 



-4'£ ( tr{P T B k P}(A k ) ii 

Jk=l 


U 


-4 X>{P T c„P}(c£ + <?„)« 

m=l 

(A.17) 

— = 4,r((S T SmP T P)„ + (PP T U + PH 

-4^2tr{S T A k S}(B k )u 

Jb=l 

-4^ir{S'-C„S}(Cl+C„) j . 

msl 

(A.18) 

The algorithm for solving (A.14) and (A.15) then takes the form 

0. ) Choose P(°> > 0, 5(°> > 0. Set k = 0. 

1. ) Set Jfc = k + 1. 

2. ) Evaluate SS^ k \ SpW using (A.16). 

3. ) Increment S W and PW as 


5(fc+i) = ,§(<«) + &$$<*) 

(A.19) 

J5(fc+1) _ p(k) + £ k $pW 

(A.20) 


where £* is chosen to ensure that each iteration brings S (fc+1) and P (fc+1) closer to convergence. 
4.) Go to 1. 

We are deliberately vague in Step 3 of the above algorithm. It would be very straightforward 
to base the choice of £ k on achieving some improvement in the cost J ; for example, choosing £ k to 
satisfy 
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J(S<* +l >,P< fc+1 >) < J(5< fc \P<*>) (A.21) 

Unfortunately, there is no information about the value of J directly available from the calculations 
in (A.16). Also, it can be seen by comparing (A.13) with (A.14 - A.18) that evaluation of (A.21) 
would add a very significant increment of numerical overhead to the algorithm, particularly when 
n is large, and n, and n p are relatively small. 

It would be much more efficient to base the test for improvement on quantities such as gradient 
norms or, better yet, the gradient traces. Work is in progress for determining an efficient, but 
reliable, improvement test. In testing the above algorithm numerically on a 20**-order problem (H 
was 400 tk -order), Step 3 was omitted entirely and £k was fixed at £ = 0.25, with good success. 

Note that the numerical overhead for evaluating (A.16) varies significantly with n, and n p . 
From (A.14 - A.16) and n, = n p = n, we see that one pass through Step 2 of the algorithm 
requires 0(4 (ij + v + 5/4)n s ) multiplications. Unfortunately, when n, < n or n p < n, the resulting 
approximation (A.2, A.12) is rank deficient. For treating cases where a sequence {Hk : k = 1,...} 
must be approximated, in which Hk is "dose” to Hk- i, an efficient algorithmic approach would 
consist of approximating Hi using n, = n p = n, then approximating subsequent members of the 
sequence for n, = n p « n by 


II kn = (/JSi + S» hS|T +1 ) ® (ph + ft+i j? +1 ) (A.22) 


S* = 3i3f + E S i S 7 (A. 23) 

i= l 

PT+Y, p i p i < a - 24 > 

3=1 

where 
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(A.25) 


Sk, Pk = arg min \\Hk - #*|| • 

S.P 

and 0 < fi < 1 is chosen to prevent the norm of Sk ® Pk from growing larger than that of Hk- 
Note that, if Si > 0 and Pi > 0, then each subsequent Sk> 0 and Pk > 0. This procedure can be 
approximated in a roundabout manner by sequentially combining the first and second algorithms in 
this Appendix, but a direct method of calculating Sk, Pk would be much more satisfying. Although 
it is unlikely that the derivation of such a method poses significant difficulty, time and the scope 
of the project did not permit pursuing this issue. 
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APPENDIX B. 


AN EFFICIENT ITERATIVE PROCEDURE FOR SOLVING 
DISCRETE LYAPUNOV EQUATIONS 

In the numerical algorithm described in Section 3, discrete Lyapunov equations (14) and (18) 
are solved for K and A k at each step in a over the integration interval [ot/,ot 0 ]. It can plainly be 
seen that, regardless of the solution procedure used, this is one of the dominant numerical overhead 
items in the algorithm. 

The most efficient procedure to date [24] for solving the equation 

A t XA-X + C = 0 (B.l) 

where C = C T is to: 

1. ) Transform A into upper real Schur form A via 

A = U t AU (B.2) 

2. ) Transform C to C via 

C = U T CU (B.3) 

3. ) Partition A according to the dimensions of the nonzero main diagonal blocks. This naturally 

breaks the problem of solving 


A t XA — X + C = 0 

into lxl and 2x2 subproblems which can be solved sequentially. 


(B.4) 
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4.) Back-transform X to X using 


X = UXU T (B.5) 

The virtue of this procedure is the operation in Step 3, since partition-wise solution of (B.4) 
requires only minimal numerical overhead, compared with solution of (B.l) as a linear system of 
order n(n + l)/2. Unfortunately, Step 1 can be quite expensive, since generally, the QR algorithm 
is used to calculate U. Because of this, the overall algorithm requires 2(3 -f 2<r)n* floating point 
multiplications, where a is the number of QR iterations required for satisfactory convergence. This 
number tends to increase when the eigenvalues of A have similar magnitudes. Because of this, when 
a solution to (B.l) is required for A which is asymptotically stable in discrete time, as in feedback 
optimization problems, a can become large. Another drawback of the above algorithm becomes 
apparent when (B.l) must be solved repetitively for a series of similar, though not identical A 
matrices. In this case, new U matrices must be calculated for each A. 

These considerations motivate investigating procedures which make use of “local” information. 
A successful algorithm of this type would exploit the fact that solution X(Ai), for A = At in (B.l), 
is “close” to X(Aj) when A\ is close to A t . A straightforward point-of-departure is to consider a 
Newton algorithm for minimizing (driving to zero) the function 

J = tr{\A T X T A -X T + Q t )[A t XA - X + Q]> (B.6) 

J in (B.6) is the squared Euclidean norm of the error in (B.l) when X is not the solution. The 
gradient of J is 

dJ/dX = 2 [AA t XAA t - A t XA - AXA t + AQA t + X - Q] (B.7) 

or 
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a J/d vec X = 2 [AA t ® AA T - A T ® A T - A® A+ I\vecX + vec [AQA T - Q] (B.8) 


and its Hessian is 


— fnF = 2[AA T ® AA T -A T ® A t -A® A+i\ (B.9) 

dvecXovec T X 1 

The Newton algorithm for this problem consists of calculating a sequence {X*} of approximations 
to the solution X f where 


vecXk+i = vecX k - (- — jrzz ) -J 

\d vec X d vec T X / dvecX I 

and 0 < £* < 1 is a parameter chosen to ensure satisfactory improvement in J from iteration to 
iteration. The algorithm, based on a Taylor series for dJ/dX about zero, does satisfy the require- 
ment for exploiting local information, but is, nonetheless, quite impractical for all but very small 
problems. This is because, for A € £ nXn , mere construction of the Hessian matrix requires 2 n 4 
multiplications, and its inversion requires 0(n®/3) multiplications. After inversion, an additional 
n 3 (n* + l)/2 multiplications are required to evaluate (B.10). The 6n* multiplications per iteration 
required to calculate dJjdX (15n s /2 on the first iteration) also impose an unattractively high 
computational burden if any significant number of iterations are anticipated. 

If we assume that solutions to (B.l) are required for a sequence of closely spaced A’s and C’s, 
a numerically inexpensive approximation to the Newton step (B.10) can be constructed, based on 
the discussion in Section A. Denote d^J/dvecXdvec^X for the i** A and C by H k in (A.22 - 
A.25). Given the form (A.22) of S k and Pk, and the fact that and are available, S'* 1 
and P can be calculated using the standard formula 


(B.10) 

x=x h 


(sw-'s’gf (s ;- 1 )*- 1 
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(SjT‘)’ = S/T-'i 


(B.12) 


where S 3 k is the j th column of Sk, and similarly for P k l . Updating S k x and P k x in this manner 
requires only 0(3n,n J ) multiplications. Since it is assumed that the A* and C* (and, hence, the 
Xk) are closely spaced, evaluation of the gradient dJ/dX cam be accomplished via numerical differ- 
encing, eliminating 6 n 3 multiplications per iteration. Overall, then, one concludes that for small n,, 
n p> only some finite multiple of n a multiplications per iteration are required. The decomposition 
giving Hk is not included in this estimate, since the derivation was not finalized, but, for small n # , 
n p , should not affect the estimate materially. 
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APPENDIX C. 


PROOFS 

This Appendix contains proofs of assertions made in Section 2 and Appendix A. 

Proof_of_ Lemma 1 : 

Consider the problem of minimizing J from (12) without any constraints on the o-dynamics 
of G. The Lagrangian for this problem is 

L= [“' e + tr{S{G,K)A^}da (C.l) 

Ja 0 

Note that, in minimizing C, G assumes the role of a control variable. The Minimum principle gives 
the necessary condition 


dM /dG = KGA.k — B t KAA k C t = 0 (C.2) 

The necessary conditions on K and A k are, again, (18) and (19). From (18, 19) and the Lemma’s 
assumptions, it is trivial to show that K and A K are fc-times continuously differentiable functions 
of a. Therefore, G*(a), given by 


G* = K + B t KAA k C t A% (C.3) 

where (.) + denotes the Moore-Penrose inverse, has a continuous k tK derivative in a. 

Denote the ij tk element of G*(a) by ff^(a), which has the Fourier series representation 


OO 

9ij{t) = a 0 + a n co3(nxt - u n ) 
n=l 


(0.4) 


with 
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t = (a — a 0 )/(ctf - a„) 


(C.5) 


over the interval [ a 0 ,ctf \ . One can form an approximation to by truncating the series (C.4) 

at the r th term: 


§i 3 (t) = a 0 + anCos{n*t - u/ n ) 


(C.6) 


so that 


»?,(<)- 5yW= 2D 


(C.7) 


n=r+l 


It is well known that the a n in (C.4) approach zero at least as rapidly as 7 /n fc+1 for some 


constant 7 -> 0; therefore, 


- 9n{t) < X) 


nssr+l 


Since k > 0, (C.8) can be rewritten 


(C.8) 


= ios> 

The truncated series (C.6) has a minimal realization as a linear system of order 2 r, so that a 
linear system of order 2 mlr suffices to realize the first r trigonometric terms in the Fourier expansion 
of each of the elements in <7*(t). Therefore, (43) holds for 




b = — ma X'Uj 

o *,j 


I 


Proofs of Lemma 2 : 

Denote a set of eigenvectors of Ai € Z pXp as v n %, n = 1, ...,p and those of B,- 
m = 1 q. Suppose that 


Bi>0 

»=i 


but that 


det(£A i } = 0 

*- 1 

Since each A,- > 0, (C.12) implies that there is some v n such that 

AiV n = 0 

for * = 1, that is, all of the Ai share a zero eigenvalue with eigenvector v n 
Ai ® Bi satisfies 


(A,‘ ® Bi)Vn ® t V m i — AiV n ® BiW m * — 0 


for m= 1, ... ,q , so that 


which contradicts 


52(A< ® Bi)v n ® 
0=1 


(C.ll), thus implying 


w mi = 0 


E * >0 

0=1 

Similar reasoning leads to 


€ as w mi , 

(c.ll) 

(C.12) 

(C.13) 

. Therefore, each 

(C.14) 
(C.15) 
(C.16) 
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i>>o 

»=x 

Note that if (C.16) holds, then for a< > 0, t = 1, t?, 


(C.17) 


^a»A» > 0 


Since P 1°) > 0 by assumption, then 


a. = tr{B.p(°>}/||P^|| a >0 | 

which, from (A.9), implies that 5 (1) > 0. This, in turn, implies that P (l * > 0, from (A.10). 
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